PRO MAKE_RADIAL_PROFILE_8UM
  
  ;;===========================================================================
  ;; Routine to create vertically-integrated radial profiles of the 8um
  ;;   radiation fields purported from the model fits made here.

  ;; Define the COMMON BLOCK
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Since we are using routines buried in BGPS_FFORE, compile now to
  ;; they are available.
  FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_MODEL2, BGPS_FFORE_SIMPLE, $
     BGPS_FFORE_SIMPLE2, BGPS_FFORE_LOOPIE, STAR_DENSITY
  RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER
  
  ;; Here's where you can change the SUBSET used (i.e. TEMP,
  ;; training set, etc.)  ALSO the SUFF is which rotation curve used
  rotc   = 0b
  subset = ''
  
  do_Tfit = 0b
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  corr = R0 / 8.5d3
  
  ;; Loop over CASES (Td=FIT)
  FOR j=0L,2 DO BEGIN
     
     dfn = './irdc_dist_model/data/case'+$
           string(j+1,format="(I0)")+$
           '_irdc_TdFIT.sav'
     restore,dfn
     
     rho_int = dblarr(!MW.NBINS*5.)
     
     ;; FOR LOOP from Z=0 to Z=large in dZ=small increments
     FOR GZ=0.d, 1000.d, 1.d DO BEGIN
        
        ;; DO vector
        R = d
        Z = d*0.d + GZ
        rho_hi = hi_density(R,Z)
        rho_h2 = h2_density(R,Z)
        rho_star = star_density(R,Z)
        
        CASE res[5] OF
           0: fj = BGPS_FFORE_MODEL(R0,res[0],res[1],res[2],res[3],rho)
           1: fj = BGPS_FFORE_SIMPLE(R0,res[0],res[1],res[2],rho)
           2: fj = BGPS_FFORE_MODEL2(R0,res[0],res[1],res[2],rho)
           3: fj = BGPS_FFORE_SIMPLE2(R0,res[0],res[1],res[2],res[3],res[6],rho)
        ENDCASE
        
        ;; To count going both up and down, 2x all but GZ = 0.
        mult = (GZ EQ 0.) ? 1.d : 2.d
        rho_int += mult * rho
        
     ENDFOR
     
     ;; plot,d/1.d3,rho_int,tit='Case #'+string(j+1,format="(I0)"),$
     ;;      xtit='Galactocentric Radius [kpc]',$
     ;;      ytit='8 um Radiation Surface Density'
     
     command = 'rad_profile_case'+string(j+1,format="(I0)")+' = rho_int'
     print,command
     err = Execute(command)
  ENDFOR
  
  Rgal = d /1.d3
  
  save,Rgal,rad_profile_case1,rad_profile_case2,rad_profile_case3,$
       filename='./irdc_dist_model/vert_int_radial_profile_8um.sav',/ver
  
END
